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non-geosynchronous orbits, and wherever a transmitter and receiver are moving relative to each other. 
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ESTIMATION OF DOPPLER SHIFT IN SEVERAL STEPS 

Field of the Invention 

This invention relates to digital filtering of nnultiple bits input signals and 
5 computation of correlation or comvolution values. 
Background of the Invention 

In formation of correlation values from digital signals, the number of 
multiplications required to form a single correlation values can be as large as 
the number P of time slots or "chips" used to define the signal, if a conventional 
10 approach is used for computation of the correlation values. In a code division 
multiple access (CDMA) approach to communications, the communication 
symbols are multiplied or "spread" using a fixed pattern of reference symbols, 
referred to as pseudo-random (PN) code, and transmitted. At the receiver end, 
an optimal solution for detecting the received symbols, and identifying the 
15 particular code used for the spreading, employs a maximum likelihood (ML) 
decision mle. The ML solution for coded symbols is a matched filter or 
correlator, with coefficients given by the same PN code, a process known as 
"despreading". Analysis of a correlation curve for the received signal and 
determination of the time shift and amplitude of its peak value allows two 
20 questions to be answered: (1) Is the reference pattern used to form the 

correlation curve likely to be the same as the spreading pattern used at the 
transmitter; and (2) If the answer to the first question is "yes", what is the most 
likely time shift associated with the received signal. 

A CDMA communication system receiver requires use of at least 3*N 
25 matched filter values, where N is the number of rake fingers used for the 

communications. TTue matched filter coefficients used for despreading account 
for the largest fraction of the gate count in a CDMA chip and have an 
associated built-in time delay for the substantial signal processing that is 
required. For example, in use of Code Division Multiple Access (CDMA) for 
30 signal recognition in a Global Positioning System (GPS), the number of chips 
used is P = 1023, which requires about 1023 adders, arranged in about ten 
layers or more, to compute a single correlation value, and this computation is 
repeated for each of the approximately 1 023 time shifts in order to compute a 
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set of values for a full correlation curve. The time delay for (a minimum of) 
ten processing layers can be many tens of nanoseconds. 

A second category where PN codes are used is in frame synchronization 
in digital communications. The transmitted infomiation symbols are divided 
5 and packed into frames that also carry header and trailer information 

(overhead), including a preamble for synchronization of the remainder of an 
asynchronously transmitted frame. A PN code or "unique word" (UW) is 
included to assist in detecting the beginning point of the payload within a frame. 
The receiver searches for the UW within a frame, and an optimal solution is an 
10 ML synchronizer. Again, this requires use of a matched filter or correlator, 
with coefficients given by the UW used in the frame. Again, the matched filter 
coefficients used for pattern recognition account for the largest fraction of the 
gate count in frame synchronization. 

What is needed is an approach that reduces the number of arithmetic 
15 operations used in formation of a correlation value used in CDMA 

communications, in frame synchronization and in any other signal processing 
procedure requiring formation of one or more correlation values. Preferably, 
the approach should have reduced complexity and should allow use of a much 
smaller gate count and have a smaller associated time delay for computation of 
20 each correlation value. Preferably, the approach should be flexible enough to 
accommodate a signal of any magnitude, using any signal representation format 
(binary, quaternary, octal, binary coded decimal, hexadecimal, etc.) and using 
any signal over-sampling rate. 
Summary of the Invention . 
25 These needs are met by the invention, which reduces the number of 

addition operations from M-N-R to a maximum of M-N, where M is the 
maximum number of orders of magnitude needed to represent a signal 
magnitude (e.g., M=l for up to 10, M=2 for up to 100, etc. in decimal base), N 
is the number of binary digits used to represent a single "order of magnitude" 
30 (e.g., N=3 for octal, N=4 for a decimal or hexadecimal order of magnitude), 
and R is the over-sampling rate used for signal reception (e^g., R = 2-10, or 
higher). The invention computes differences between two consecutive 
correlation values, rather than computing each correlation value separately, and 
this allows reduction, by a multiplicative factor of at least R, depending upon 
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the structure of the reference or coding word, in the number of arithmetic 
operations for computation of a correlation value. The number of gates 
required is reduced by the same factor, and the associated time delay for 
computation of a correlation value is reduced by at least At(add)-log2(R), 

5 where At(add) is the time required to perform an addition of two values. 
Optionally, a time shift point at which the correlation curve reaches a 
maximum is automatically identified by this approach, and the corresponding 
(maximum) correlation value is quickly recovered after this time shift point has 
been identified. 

10 Brief Description of the Drawings 

Figure 1 is a graphical view of a representative correlation curve for a 
received signal. 

Figure 2 illustrates a gate array used to compute correlation values 
according to a conventional approach. 
15 Figure 3 illustrates overlapping arrays of over-sampled PN Code values 

(±1) used in an example that explains use of the invention. . 

Figure 4 illustrates a gate array used to compute correlation values 
according to the invention, corresponding to the conventional approach in 
Figure 2. 

20 Detailed Description of the Invention 

Figure 1 is a graph of an idealized autocorrelation function that might be 
formed from a received digital signal s(t;rec), which is generally expressible as 
a multi-bit numerical value. Each time point t^ corresponds to a time shift for a 
convolution value C(tn), computed as 

25 K 

C(tn) = Z s(tk;rec)-s(tn- tk;ref), (1) 

k=l 

where s(t;ref) is the reference signal used at the receiver to implement pattern 
recognition for the received signal and {ty^} is a sequence of time points (not 

30 necessarily equidistantly spaced) used to compute a correlation value for the 

received signal. The time shift At(max) corresponding to the position where the 
maximum amplitude of the correlation curve in Figure 1 may be, but need not 
be, precisely equal to one of the time shift values t^ used to compute the 

convolution values. 
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Use of the invention is first illustrated by an example. A PN code signal 
has an amplitude s(t;ref) = 789, which is represented in hexadecimal format as 
an ordered sequence 01 11/1000/1001. The received signal is over-sampled at a 
selected rate R, for example, R = 4, so that each of the 12 signal values 

5 representing the received signal s(t;rec) is appears four consecutive times. A 
"modified hexadecimal format", in which the binary value "0" is replaced by 
the numerical value "-1", is adopted here to more easily illustrate operation of 
the invention. A modified form of the invention will also work with the 
conventional hexadecimal format representation. Adopting the modified 

10 hexadecimal format, the over-sampled expression for PN Code amplitude 
s(t;ref) == 789 becomes 

{sCt;ref)} =-1-1-1-11111111111111111-1-1-1-1 

-1 -1 -1 -1 -1 -1 -1-11111 -1 -1 -1 -1 -1 -1-1-1111 1. (2) 
This number has 48 = 3-4-4 binary digits in its representation, including factors 

15 corresponding to an order of magnitude number, M = 3, a hexadecimal (or 

decimal) representation number, N = 4, and an over-sampling rate number, R = 
4, in this example. The time shifts tj^ are equidistant so that tj^+j - 1 = At = 

constant.. 

In a conventional approach to computation of the correlation values using 
20 this PN Code, a sequence of time shifts tj^ is applied to the over-sampled 

reference signal values s(ti^;ref), and the resulting signal values are multiplied 

by the received signal values, as indicated in Eq. (1). For any selected value of 
the time shift tj^, computation of the corresponding correlation value requires 

use of 48 signal value multipliers and use of 47 (or 48) adders, which may be 
25 arranged six layers containing 24, 12, 6, 3, 2 and 1 two-input adders as shown 
in Figure 2. The associated time delay At(conv) for conventionally processing a 
received signal to provide a single correlation value is six times the time delay 
At(add) associated with addition of two values at a single adder. More 
generally, the associated time delay is approximately {log2[M-N-R])+, where 

30 {P)^ is the smallest integer that is greater than or equal to the real number P; 
this is related to, but not the same as, the integer part, [P]^., of the real number 
P. For this example, the associated time delay is At(add)-{log2[3-4-4] }+, = 
At(add)-{ 5.585 }+ = 6-At(add). 
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For the example chosen, the convolution values C(i^) defined in Eq. (1) 
become 

K 

C(n-At) = X s(k-At;rec)-PN((n-k)-At), (3) 
5 k=l 

where the sequence of PN code values is set forth in Eq. (2). The difference of 
two consecutive convolution values is verified to be 

AC(tn+i) = C(tn+i)-C(tn) 

47 

10 =1 s(k-At;rec)-{PN((n+l-k)-At) - PN((n-k)-At) } 

k=l 

+ s(48-At;rec)-PN((n+l-48)-At) - s(0-At;rec)-PN(n-At) 

48 

= X s(k-At;rec)-w(n;k), (4) 
15 k=0 

w(n;k) = 0 (n-k not divisible by 4) 
= PN((n+l-48)-At) (k = 48) 
= -PN(n-At) (k = 0) 

= 2 (n-k = 4-u; PN(4-u) = -PN(4-u + 1) = -1; k 0, 48) 
20 =0 (n-k = 4-u; PN(4-u) = PN(4-u -H 1) ) 

= -2 (n-k = 4-u; PN(4-u) = -PN(4-u + 1) = +1), (5) 
where u = 0, ±1, ±2, ±3, ... . One verifies that the convolution value difference, 
C((n-hl)-At) - C(n-At), includes at most 13 = 3-4 + 1 terms, corresponding to 
the integer values of k in the convolution sum in Eq. (3) for which the integer 
25 difference, n-k, is divisible by 4 (n-k = 0, 1, 2, 48). In the particular 

example with a PN code of 789, of the indices for which n-k is divisible by 4, a 
sequence of 1 3 consecutive values of the coefficients w(n;k) becomes 

{w(n;k)} = {-1, 2,0,0, -2,0,0,0, 2, -2,0, 2, 1} (6) 
with the first and thirteenth values w(n;k) (= -i-l or -1) being determined by 
30 individual PN Code values, rather than Code differences. 
This example has used the particular choices 
M = maximum order of magnitude of reference signal = 3, 
N = number of binary digits used for each order of magnitude = 4, 
R = over-sampling rate = 4. 
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Figure 3 illustrates the 48 values of the PN Code represented in Eq. (2), 
for the shifted time values s(t;ref), s(t+At;ref) and s(t+2At;ref). Each of these 
three sequences of 48 numbers is multiplied by corresponding members of a 
digital signal sequence s(t;rec), as indicated in Eq. (1). Because of the over- 

5 sampling rate (R=4) adopted in this example, every fourth value of the 
sequence in Figure 3 representing s(t;ref) is a potential "boundary point", 
where the PN Code sequence members s(t;ref) and s(t+At;ref) can have 
different signs (e.g., +1/-1 or -1/+1)- Similarly, every fourth value of the 
sequence representing s(t+At;ref) is a potential "boundary point", where 

10 s(t+At;ref) and s(t+2At;ref) have different signs. Where two time-shifted 

sequences, such as s(t;ref) and s(t+At;ref) have the same sign (e.g., -1-I/-1-I or 
-1/-1) at such a boundary point, this boundary point will not contribute a non- 
zero value to a convolution value difference C(tj^) - C(tj^.i). 

Figure 4 illustrates an array 41 of coefficient multipliers and signal 
15 summers that can be used to form a convolution value difference for the 

example considered here, with M = 3, N = 4 and R = 4 but with an arbitrary 
PN code. The top level of the array has a maximum of M-N-1 = 1 1 coefficient 
multiplier modules 43-k, each with an associated coefficient w(n;k) for k = 1, 

2 11 that has a value -i-2, 0, or -2, depending upon the particular PN Code 

20 used, plus one signal differencer that receives the signals s(48-At;rec) and 
s(0;rec) and forms and issues a weighted difference 

As(n) = s(48-At;rec)-PN((n+l-48)-At) - s(0;rec)-PN(n-At). (7) 
The output signal for the coefficient multipher module 43-k is 

0(43-k) = w(n;k)-s(4-k-At;rec). (8) 
25 For the PN Code 789 used in the preceding example, 
w(n;l) = w(n;8) = w(n;ll) = +2, 
w(n;4) = w(n;10) = -2. 

w(n;2) = w(n;3) = w(n;5) = w(n;6) = w(n;7) = w(n;9) = 0, 
as can be verified from Eq. (6). The second array level from the top has six 
30 two-input signal adders 45-j (j = 1, 2, 6) producing six output sum signals. 
The third array level from the top has three two-input signal adders 47-h (h = 
1, 2, 3), producing three output sum signals. The fourth array level from the 
top has one two-input adder 49; and the fifth array level from the top has one 
two-input signal adder 51, whose output signal is the convolution value 
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difference AC(ti-,) for some index value n. The array 41 in Figure 4 has M-N-1 
= 1 1 signal coefficient multipliers at the top level and four levels with an 
additional 12 signal adders. The maximum associated time delay is 
approximately 5 At(add). The number of gates required is a maximum of 23 
5 (18 for the PN Code 789 chosen for this example), as compared with 47 for the 
conventional approach for convolution. 

More generally, the convolution value set forth in Eq. (4) can be 

expressed as 

AC(tn+i) = C(tn+i)-C(tn) 

10 M-N-R 

= X s(k-At;rec)-{PN((n+l-k)-At) - PN((n-k)-At)} 

k=0 
M-N-R 

= Z s(k-At;rec)-W(n;k), (9) 
15 k=0, 

W(n;k) = 0 ( n-k not divisible by R) 

= PN((n+l -M-N-R)- At) (k = M-N-R) 
= -PN(n-At) (k = 0) 

= 2 ( n-k = R-u; PN(R-u) = -PN(R-u -t- 1) = -1) 
20 =0 (n-k = R-u; PN(R-u) = PN(R-u+ 1)) 

= -2 ( n-k = R-u; PN(R-u) = -PN(R-u + 1) = +1), (10) 
and the maximum number of terms to be added to compute the convolution 
difference is M-N + 1 . MultipHcation by a factor of 2 is implemented using a 
one -bit shift of a binary-expressed value. 

The innovative procedure begins with computation or provision of a first 
convolution value 

M-N-R 

C(ti) = L s(tk;rec)-s(ti- tk;ref). (1 1) 

k=0 

.^0 This first convolution value may be known from other considerations (e.g., 
C(ti) = 0) and may not need to be computed. For n > 1, only the convolution 
differences, AC(tn+i) = C(tn+i) - C(tn), as set forth in Eq. (9), are computed. 
The total number of additions (and multiplications by 2) required here is at 
most M-N and may be less, depending upon the values of the coefficients 



25 
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W(n;k) as set forth in Eq. (10). One can verify that the number of layers of 
two-input adders needed to accomplish these M-N value additions is at most 
{log2[M-N}^. If K-input adders (K>2) are substituted for the two-input adders, 

the number of layers of K-input adders needed to accomplish these M*N value 
5 additions is reduced to at most {logj^[M-N}4.. Each factor-of-two 

multiplication required by the coefficients W(n;k) in Eq. (10) can be 
implemented by shifting left by one binary digit in a shift register. If a 
particular convolution value, say C(tj^), is needed, this value is computed using 

the relation 
10 m 

C(tm) = C(ti) + X {C(tn) - C(tn.i)). (12) 

n=l 

For the particular example discussed in the preceding. Figure 4 illustrates 
a suitable array to compute the convolution value differences. This array 
15 requires at most {log2[M-N}+ = 4 layers, as opposed to {log2[M-N-R]}^, = 6 

for the conventional approach, and requires at most M-N =12 adders, as 
opposed to M-N-R = 48 for the conventional approach. 

The integer M depends upon the "order of magnitude" chosen for the PN 
code. The preceding discussion relies upon use of a hexadecimal format, in 

20 which N=4 consecutive binary digits are used to express each order of 

magnitude value, such as a decimal. Use of a decimal format would also require 
N=4, because use of at least four binary digits is required to express the values 
8 and 9 in a binary system. One can also express an order of magnitude value in 
octal format, using N=3 consecutive binary digits. In an octal format, the 

25 largest order of magnitude value becomes 7 (expressed as 111 in an octal 

format), not the maximum order of magnitude value 9 of a decimal format or 

the value F (equivalent to 15) of a hexadecimal format. 

The innovative method disclosed here uses single unit differences 
AC(tj^^i) = C(tj^^.j) - C(tj^), as set forth in Eq. (9), to quickly compute the 

30 convolution values themselves, as set forth in Eq. (12). Double unit differences 

^2C(tn+2) = C(tn+2) ' C(tn) 

= AC(tn+2) + ^C(tn+i) (13) 
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can be computed in a similar manner, if desired, using a sum of single unit 
differences, as computed in the preceding. More generally, one can compute p- 
unit differences 

ApC(tn+p) = C(tn+2) - C(tn) 

5 P 

= ZAC(tn+k) (14) 
k=l 

in a similar manner, if desired. 

The innovative approach disclosed here is also useful in estimating the 
10 time shift value At(max) that corresponds to a local maximum value of the 
convolution value C(t,-,), where At(max) need not coincide with one of the 
discrete time shift values tp but may lie between two discrete time shift values. 
Assume that tji^.j , ij^, and t^^+i , are three consecutive time shift points for the 

correlation curve in Figure 1 , and that the time shift value At(max) lies between 
^5 tm-1 ^m+1- situation assumes that the consecutive convolution values 
C(tin_i), C(tm) and C(tm+i) satisfy the constraint 

C(tn,) >max{C(tm-i), C(tm+i)}- (15) 
This constraint may be expressed in an equivalent form for convolution signal 

differences as a pair of constraints 
20 AC(tm)^0, (16A) 

AC(tm-M)^0. (16B) 

The value T(max) can be computed as follows. A second degree (or higher 
degree) polynomial fit to the correlation values at the time shift points tfj^.j, 

tm, and tj^+i , becomes 
25 P(t;2) = cO + c1'T + c2-x2 (tj^.j <T<tm+l), (17) 

c2 = {lC(tm.i) - C(tm+i)]-(tm - tm+l) 

+ lC(tm) - C(tin+i)]-(tm+l - tm-l)}/(det), (18) 

cl = {[C(tm-i) - C(tm+i)]-(tn,+ l2- tm^) 

+ [C(tm) - C(tn,+i)]-(tm.i2- tm+i2)}/(det), (19) 

30 cO = C(tm)-cl-(tm)-c2-(tm)2, (20) 

det = {(tm.i -tmHtm -tm+lHWl "Wl))"^- (21) 

The time shift for maximum amplitude is given by 

T(max) = -cl/c2, (22) 
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and the time shift value T(max) depends only upon the two correlation signal 
differences, AC(tm) = C(t^) - C(Atm-i) - C(Atm+l) and AC(tni+i) = 
C(Atj^4.l) - C(Atjn), which have already been computed according to the 

procedure set forth in the preceding, 
5 More generally, the second degree polynomial P(x;2) can be replaced by 

an nth degree polynomial P(x;n) (n>2), with appropriate constraints imposed on 
the polynomial values at three or more of the discrete time points {tj. }, and the 

time shift value T(max) that maximizes the polynomial P(T;n) within the 
interval tj^^.j < x < tj^+l can be determined. 

10 In one interpretation, the time shift value x(max) is identified as an 

estimate of the time shift value corresponding to the maximum value the 
convolution value C(tj^) would attain, if the convolution were extended to 

continuous values of the time variable t. 

In a second interpretation, one imposes the conditions (16A) and (16B) 
15 and determines or computes a time variable value t = tm' for which 

Itm* - = mini It^' - tm-il, Itm' - tml, Itm' " Wll >• (^3) 
The time variable value tj^» is then interpreted as the time shift value, drawn 
from the discrete set {tj^}, that produces the maximum convolution value. 

In a third interpretation, one imposes the conditions (16A) and (16B) and 
20 determines or computes a time variable value t = tj^" for which 

C(tm") = max{ C(tm-l), C(tm), Ca^+i) 1, (24) 
which will normally be t^^" = tj^^. 

The time shift value t = x(max) or t = tj^^ or t = tj^" can be used to 

determine the most likely time shift in CDMA communications, the most likely 
25 point at which the payload (data) begins in frame synchronization, or for any 
other similar purpose involving use of convolution signals. 
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1. A method for carrier signal recovery in the presence of a Doppler 
shift, the method comprising: 

providing a carrier signal that may contain a Doppler-shifted replica of a 
desired signal; 

processing a first portion of the carrier signal through a first linear 
predictor to obtain a first contribution to Doppler frequency offset for the 
received signal; 

modifying a second portion of the carrier signal to obtain a first- 
modified signal from which the first offset contribution is removed; 

processing the first-modified signal through a second linear predictor to 
obtain a second contribution to Doppler frequency offset for the received 
signal; 

modifying the first-modified signal to obtain a second-modified signal 
from which the second offset contribution is removed; and 

processing the second-modified signal through a phase locked loop circuit 
to estimate at least one of a residual contribution to Doppler-frequency offset 
and a phase angle for the provided signal. 

2. The method of claim 1 , further comprising estimating said Doppler 
frequency offset for said carrier signal as a sum of said first, second and third 
Doppler frequency offset contributions. 

3. The method of claim 1, further comprising providing as said first 
portion of said carrier signal a carrier signal segment including at most 10 
symbols. 

4. The method of claim 3, further comprising providing as said second 
portion of said carrier signal a second carrier signal segment including at most 
30 symbols. 

5. The method of claim 1, wherein said process of modifying said second 
portion of said carrier signal to obtain said first-modified signal comprises 
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multiplying said second portion of said carrier signal by a sinusoidal signal 
having a frequency equal to said first contribution to said Doppler frequency 
offset. 



6. The method of claim 5, wherein said process of modifying said first- 
modified signal to obtain said second-modified signal comprises multiplying 
said first-modified signal by a sinusoidal signal having a frequency equal to said 
second contribution to said Doppler frequency offset. 

7. The method of claim 1, wherein said step of processing said second 
modified signal through said phase locked loop comprises processing said 
second modified signal through a decision feedback phase locked loop to 
estimate at least one of said residual contribution to said Doppler-frequency 
offset and said phase angle for said provided signal. 

8. A method for carrier signal recovery in the presence of a Doppler 
shift, the method comprising: 

providing a first carrier signal sl(t) that may contain a Doppler-shifted 
replica of a desired signal; 

computing a first correlation function 

Nl 

Sl(Nl;Kl;Ts) = Z sl(k-Tjj;carr)-sl((Kl-k>Ts;carr)* 
k=l 

of the carrier signal with itself, where Nl (>1) is a first selected sample size, 
Kl-Tg is a first selected time delay and I/T3 is a selected signal sampling rate; 

computing a first Doppler frequency offset component 
fDi =tan-l{Sl(Nl;Kl;Ts))/(27i-Kl-Ts); 

forming a first-modified carrier signal 

s2(t;carr) = sl(t;carr) exp(-j27C-f£)-j*t); 

computing a second correlation function 

N2 

S2(N2;K2;Ts) = Z s2(k-Ts;carr)-s2((k-K2)-Ts;carr)* 
k=l 
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of the modified carrier signal with itself, where N2 (>1) is a second selected 
sample size and Kl-T^ is a second selected time delay; 

computing a second Doppler frequency offset component 
fD2 = tan-1 {S2(N2;K2;Ts)}/(27r-K2-Ts); 

forming a second-modified carrier signal 

s3(t;carr) = s2(t;carr) exp("j27i;-fj)2*t); and 

processing the second-modified signal through a phase locked loop circuit 
to estimate at least one of a residual contribution to a Doppler-frequency offset 
and phase angle for the provided signal. 

9. The method of claim 8, further comprising estimating said Doppler 
frequency offset for said received carrier signal as a sum of said first, second 
and third Doppler frequency offsets. 

10* The method of claim 8, further comprising: 
forming a second modified carrier signal 

s3(t;carr) = s2(t;carr) exp(-j27C-fD2*t)r and 

processing the second modified carrier signal by feedback decision phase 
locked loop analysis to obtain an estimate of at least one of a third Doppler 
frequency shift and a phase associated with the second modified carrier signal. 

1 I. Tlie method of claim 8, further comprising choosing said integers Nl 
and N2 to satisfy the relation Nl < N2. 

12. The method of claim 1 1, further comprising choosing said integers 
Nl and N2 to be Nl = 60 and N2 = 180, 

13. The method of claim 8, further comprising choosing said numbers 
Kl and K2 to satisfy the relation Dl < D2. 

14. The method of claim 8, further comprising choosing said numbers 
Kl and K2 to be Kl = 1 and K2 = 20. 
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15. The method of claim 8, wherein said step of processing said second- 
modified signal through said phase locked loop circuit further comprises 
processing said second-modified signal through a decision feedback phase 
locked loop circuit to estimate at least one of said residual contribution to said 
Doppler frequency offset and said phase angle for said provided signal. 

16. A system for carrier signal recovery in the presence of a Doppler 
shift, the system comprising: 

a signal source that provides a carrier signal that may contain a Doppler- 
shifted replica of a desired signal; and 
a computer that is programmed: 



predictor to obtain a first contribution to Doppler frequency offset for the 
provided signal; 



modified signal from which the first offset contribution is removed; 

to process the first-modified signal through a second linear 
predictor to obtain a second contribution to Doppler frequency offset for the 
provided signal; 

to modify the first-modified signal to obtain a second-modified 
signal from which the second offset contribution is removed; and 

to process the second-modified signal through a phase locked loop 
circuit to estimate at least one of a residual contribution to Doppler-frequency 
offset and a phase angle for the provided signal. 

17. The system of claim 16, wherein said computer is further 
programmed to estimate said Doppler frequency offset for said carrier signal as 
a sum of said first, second and third Doppler frequency offset contributions. 

18. The system of claim 16, wherein said computer is further 
programmed to provide as said first portion of said carrier signal a carrier 
signal segment including at most 10 symbols. 



to process a first portion of the carrier signal through a first linear 



to modify a second portion of the carrier signal to obtain a first- 
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19. The system of claim 18, wherein said computer is further 
programmed to provide as said second portion of said carrier signal a second 
carrier signal segment including at most 30 symbols. 

20. The system of claim 16, wherein said computer is further 
programmed to modify said second portion of said carrier signal to obtain said 
first-modified signal by multiplying said second portion of said carrier signal 
by a sinusoidal signal having a frequency equal to said first contribution to said 
Doppler frequency offset. 

21. The method of claim 20, wherein said computer is further 
programmed to modify said first-modified signal to obtain said second- 
modified signal by multiplying said first-modified signal by a sinusoidal signal 
having a frequency equal to said second contribution to said Doppler frequency 
offset. 



22, The system of claim 16, wherein said computer is further 
programmed to process said second modified signal through said phase locked 
loop by processing said second modified signal through a decision feedback 
phase locked loop to estimate at least one of said residual contribution to said 
Doppler-frequency offset and said phase angle for said provided signal. 

23. The system of claim 16, wherein said computer is further 
prograrmned to process said first portion of the carrier signal through a first 
linear predictor to obtain a first contribution to Doppler frequency offset for 
the received signal by: 

computing a first correlation function 

Nl 

S1(N1;K1;T3) = I sl(k-Ts;carr>sl((Kl-k)-Ts;carr)* 
k=:l 

of the provided with itself, where Nl (>1) is a first selected sample size, Kl-T^ 
is a first selected time delay and l/T^ is a selected signal sampling rate; and 

computing a first Doppler frequency offset component 
f D 1 = tan- 1 { S 1 (N 1 ;K 1 ;Ts) } /(27r-K 1 -T^). 
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24. The system of claim 23, wherein said computer is further 
programmed to process said second portion of said provided signal through a 
second linear predictor to obtain a second contribution to Doppler frequency 
offset for the received signal by: 

forming a first-modified signal 

s2(t;carr) = sl(t;carr) exp("j27C-fDi -t); 

computing a second correlation function 

N2 

S2(N2;K2;Ts) = L s2(k-Ts;carr)-s2((k-K2)-Ts;carr)* 
k==l 

of said first-modified signal with itself, where N2 (>1) is a second selected 
sample size and K2-Tg is a second selected time delay; and 

computing a second Doppler frequency offset component 
fD2 = tan-MS2(N2;K2;Ts)}/(27C-K2-Ts). 

25. The system of claim 24, further comprising forming said second- 
modified carrier signal s3(t;carr) according to 

s3(t;carr) = s2(t;carr) exp(-j27C-fD2'0- 
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